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We investigate the evolution of a light impurity particle in a Lorentz gas where the background 
atoms are in thermal equilibrium. As in the standard Lorentz gas, we assume that the particle is 
negligibly light in comparison with the background atoms. The thermal motion of atoms causes 
the average particle speed to grow. In the case of the hard-sphere particle-atom interaction, the 
temporal growth is ballistic, while generally it is sub-linear. For the particle-atom potential that 
diverges as r~ x in the small separation limit, the average particle speed grows as *Vf 2 (d-i)+A) in d 
dimensions. The particle displacement exhibits a universal growth, linear in time and the average 
(thermal) speed of the atoms. Surprisingly, the asymptotic growth is independent of the gas density 
and the particle-atom interaction. The velocity and position distributions approach universal scaling 
forms which are non-Gaussian. We determine the velocity distribution in arbitrary dimension and 
for arbitrary interaction exponent A. For the hard-sphere particle-atom interaction, we compute the 
position distribution and the joint velocity-position distribution. 

PACS numbers: 05.20.Dd: Kinetic theory, 45.50.Tn: Collisions, 05.60.-k: Transport processes 



I. INTRODUCTION 

The goal of this work is to investigate the behavior of 
an impurity particle (particle in short) in a monoatomic 
gas. We focus on the limit when the particle is negligibly 
light in comparison with background atoms. In other 
words, the particle is affected by collisions with atoms, 
while atoms do not "feel" the presence of the particle. We 
want to understand the evolution of the particle velocity 
and displacement distribution. 

The problem is a natural generalization of the stan- 
dard Lorentz gas [THS] where scatters are assumed to be 
immobile. The speed of the particle remains constant 
in the framework of the Lorentz model. In our model 
the behavior is completely different and can be simply 
understood using arguments from the equipartition the- 
orem (when the background gas has a positive tempera- 
ture the average speed of the particle increases without a 
bound since the particle "tries" to reach an equilibrium 
with the background atoms). 

The problem is also reminiscent of the model originally 
proposed by Fermi [7], and later refined by Ulam [8], 
to explain the acceleration of interstellar particles and 
cosmic rays. Fermi's acceleration mechanism has been 
mostly studied using methods of dynamical systems (see 
[9] and references therein) ; an application of kinetic the- 
ory to Fermi's mechanism has been presented in |10) . 

Here we analyze the behavior of the light particle in 
an equilibrium gas using the Boltzmann equation frame- 
work. The Boltzmann equation is the basic tool in 
elucidating the properties of transport phenomena. The 
non-linear integro-differential Boltzmann equation is so 
formidable, however, that apart from the equilibrium 
Maxwell-Boltzmann distribution |T2] there are essentially 
no solutions to the Boltzmann equation [13j ■ The stan- 
dard Lorentz gas model where a point particle is elas- 
tically scattered by immobile hard spheres is described 
by the Lorentz-Boltzmann equation [1, which is linear 



and, not surprisingly, amenable to analytical treatments. 
The Lorentz gas has played an outstanding role in con- 
crete calculations (e.g. of the diffusion coefficient) and 
in the conceptual development of kinetic theory [21 13]. 
Yet the very applicability of the Boltzmann framework 
to the Lorentz gas is questionable — when the scatters 
are fixed, the molecular chaos assumption underlying the 
Boltzmann equation cannot be justified [2f(6]. 

If, however, the background atoms move and collide 
with each other, the molecular chaos assumption holds in 
the dilute limit and the (properly generalized) Lorentz- 
Boltzmann equation must be applicable as long as the 
mass of the particle is infinitesimally small so that it 
does not affect the motion of atoms. Moreover, since 
the (average) particle speed continues to grow, it eventu- 
ally greatly exceeds the typical velocities of background 
atoms. This allows to simplify the most difficult term in 
the Boltzmann equation, the so-called collision integral; 
mathematically, an integral operator becomes a differen- 
tial one and the integro-differential Lorentz-Boltzmann 
equation reduces to a partial differential equation. 

The unlimited velocity growth suggests that the parti- 
cle velocity distribution approaches a scaling form. The 
scaled velocity distribution satisfies an ordinary differen- 
tial equation (Sects. |ll}|IV| ) which admits a simple solu- 
tion; for the hard-sphere atoms, the scaled velocity dis- 
tribution is exponential (Sects. [Il III I. The Boltzmann 
equation approach also describes the spatial distribution 
of the particle, yet extracting the density distribution is 
much more difficult as it does not obey a closed equa- 
tion, so one must rely on the joint distribution function 
that simultaneously describes the probability density for 
the position and velocity. In Sec. [V] we outline the evo- 
lution of the displacement using heuristic arguments and 
exact calculations in one dimension based on the veloc- 
ity correlation functions. In Sect. VI wc derive kinetic 



equations describing the joint distribution in the long- 
time limit. In Sect. |Vli| we investigate the density profile 
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of the hard-sphere gas by utilizing the moment approach 
and in Sect. VIII we compute the joint distribution. We 
report the results of numerical simulations in Sec. |IX| and 
summarize our findings in Sect. |X| 



II. ONE DIMENSION 

As a warm-up, consider the one-dimensional case. This 
may appear physically dubious as the particle is caged 
between two adjacent atoms, so the molecular chaos as- 
sumption (that is, the lack of correlations between pre- 
collision velocities) underlying the Boltzmann approach 
is certainly invalid in one dimension. A Boltzmann equa- 
tion, however, makes sense if we consider the situation 
when in each collision the scattering occurs with a cer- 
tain probability (otherwise the particle and an atom just 
pass through each other). This one-dimensional Boltz- 
mann equation sheds light on the three-dimensional case. 
Therefore it has been proven useful as a toy model and it 
has been studied in a number of one-dimensional settings 
(see e.g. PHM3). 

The Boltzmann equation for the particle velocity dis- 
tribution f(v,t) reads 



df(v,t) 
dt 



du\v-u\P(u)[f(2u-v,t)-f(v,t)] (1) 



Here 



P(u) = 



2 /2T 



(2) 



is the equilibrium velocity distribution of the background 
atoms corresponding to temperature T (we set the atomic 
mass to unity). We shall see, however, that we do not 
need the detailed form ^ of the equilibrium Maxwell- 
Boltzmann distribution. To establish the asymptotic be- 
havior of f(v, t) it is sufficient to assume that P{u) is an 
even function, P(u) = P(—u). Even a weaker condition 
that the average velocity of atoms vanishes, 



duuP{u) = 0, 



(3) 



suffices. Whenever ^ holds, the long-time behavior de- 
pends only on the second moment of P{u) which essen- 
tially defines the temperature: 



du u 2 P(u) = pT 



(4) 



We shall see that in the long-time, more precisely when 

t » p^T- 1 ' 2 (5) 

the Boltzmann equation for the particle velocity dis- 
tribution simplifies to 



df df , d 2 f 



dv 



dv 2 



2pTt 



(6) 



This kinetic equation admits the scaling solution 



f(v,t) 



1 

2r 



-\v\h 



(7) 



To derive (Js])— ([7]) we first simplify the collision integral 
in Eq. ([lj in the t — > oo limit. Since f(v : t) = f(—v, t), it 
suffices to investigate the v > region [TS] . Moreover we 
can replace \v — u\ by v — u since the region v < u where 
the replacement is invalid provides a negligible contribu- 
tion in the long-time limit: P(u) is very small in this 
region. More precisely, the above simplification applies 
if the average speed of atoms (it) ~ VT is much smaller 
than the particle velocity v. This is our working assump- 
tion which will be checked a posteriori. When (u) <C v 
we can additionally expand f(2u—v) that appears in the 
collision integral in Eq. ([I]) into a Taylor series 



f(2 U -v)=f(v)-2u 



df(v) 
dv 



2u 



2 d 2 f(v) 



dv 2 



(2u) 3 d 3 f(v) + (2m) 4 8 4 f(v) 



3! 



dv 3 



4! dv 4 



Plugging this expansion into Eq. ([!]) and computing the 
integrals over u we obtain 



df 
dr 



df 

dv 



d 2 f fd 3 f d A f 

— - -+- 2T I — - -+- 7. — J - 



dv 2 



dv 4 



(8) 



In computing the integrals leading to the first two terms 
on the right-hand side of Q it suffices to use the integral 
relations ([3])-(|4]). The next two terms are obtained using 
the integral relations 



duu 3 P(u) = 0, 



du u A P(u) = ZpT 2 (9) 



The first relation in (pM is valid for any symmetric ve- 
locity distribution, P(u) = P(—u), while the second is 
derived from the equilibrium Maxwell-Boltzmann distri- 
bution p). 

The first two terms on the right-hand side of ^ scale 
as t _1 , the next two terms scale as Tt~ 3 , so they are 
asymptotically negligible when r 3> VT, that is, the av- 
erage particle speed greatly exceeds the average speed of 
atoms. The two following terms [which haven't been dis- 
played in Q] contains T 2 ^f and T 2 u^^, so they scale 
as T 2 t~ 5 and therefore they are even smaller. Thus in 
the t 3> VT limit (which is given by Eq. ^ in the orig- 
inal variables), Eq. Q indeed reduces to Eq. ^ in the 
leading order. 

The form of equation ^ suggests to seek the scaling 
solution of the form 



f(v,r) 



w 



v/t 



(10) 



Plugging ( 10 1 into (JgJ) we obtain an ordinary differential 
equation for $(iu) which is solved to yield $>(w) — C e~ w . 
Recalling that the particle velocity distribution is even 
and using the normalization condition f dv f{v,t) = 1 



3 



fixes the amplitude C = 1/2 and leads to the announced 
result Q. 

Having determined the scaling solution j7j|, we would 
like to understand if any arbitrary function f(v,t) ap- 
proach the scaling solution |7| in the long time limit. 
The answer to this question is presumably affirmative, 
at least when the initial velocity distribution f(v,t = 0) 
quickly decays when \v\ —¥ oo. Yet to prove this asser- 
tion even for simplest initial velocity distributions like 
f(v,t = 0) = S(v) is hard. Analytical arguments show- 
ing that the scaling solution Q is indeed an attractor 
are presented in Appendix |A") 



III. HARD-SPHERE GAS 

Consider now the most natural three-dimensional sit- 
uation and assume that atoms are hard spheres of radius 
a. We ignore both the mass and the size of the particle. 
The latter assumption is not crucial — if the particle is 
a sphere of radius b, it suffices to replace a by a + b in 
the following formulae. 

We again employ the Boltzmann equation approach. 
This framework is applicable only in the diluted limit; for 
the hard-sphere gas, this means that the volume fraction 
occupied by atoms is small: p x y« 3 < 1 (here p is the 
number density of background atoms). 

The Boltzmann equation reads 



a/(y,*) 

dt 



= j 'duP(u)ga 2 Jl)e[f(V,t)- f(y,t)] (11) 



Here e is the unit vector pointing to the position of the 
particle at the moment when it hits the sphere. The 
post-collision velocity v' of the particle can be expressed 
via v, e, and the relative velocity g = u — v: 



v' = v + 2e(g • e) 



(12) 



In Eq. (Ill we have also used the shorthand notation De 



for the integration measure over angular coordinates. For 
the hard-sphere gas, this integration measure reads [3] 



(gig) 
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(g • e) d 2 e 



(13) 



In the above expression #(•) is the Heaviside step function 
and d 2 e is the standard angular integration measure. 

we shall pro- 



To simplify the Boltzmann equation ( 11 



ceed as in one dimension. Since the partic e velocity dis- 
tribution is (asymptotically) isotropic, let us treat /(v) 
as a function of V = v 2 = (v • v). Squaring |l2]) we get 



V = V + 4(v e)(g • e) + 4(g • e) 2 = V + 4(u • e)(g • e) 

Using this result and expanding /(v') = f(V) into a 
Taylor series we obtain 



f(V) = /(U)+4(u-e)(g.e) |£ + 8(u-e) 2 (g.e 



d 2 / 
dV 2 



Using this expansion we simplify (111 to 



dl_ A df 

dt 



4 wJ duP{u) ga2 J De (u ' e) (g ' e) 

8^JduP(u)ga 2 J De (u • e) 2 (g • e) 



(14) 



As in the one-dimensional case, it suffices to keep only 
the terms with the first and second order derivatives in 
V; the terms with higher order derivatives are asymptot- 
ically negligible. The angular integrals in Eq. (14) are 
computed [see Appendix [B] to yield 



J Be(u • e)(g • e) = ~ (g • u) 
J 2e(u.e) 2 (g.e) 2 ^ [3(g-u) 



2 i 2 21 

+ g u \ 



(15a) 
(15b) 



Inserting (15a|-(15b) into Eq. (14 1 we obtain 



1 df v d 2 f 



2na 2 dt 3 dV 2 
dV 



J duP(u) [3(v • u) 2 +v 2 u 2 ] 
JduP(u)g(g ■ u) 



(16) 



In the first integral we already replaced g by — v which 
is correct in the leading order. In the second integral we 
should be more careful. We write 

g(g ■ u) = -v(v ■ u) + v^ 1 [(v • u) 2 + v 2 u 2 ] + . . . 

The integral that contains the leading term vanishes since 
/ rfuf(u)u = 0. Thus Eq. (filj} becomes 



1 df v d 2 f 



2na 2 dt 3 dV 2 



J duP(u) [3(v-u) 2 +i;V] 



(17) 



Using relations 
J e?uP(u) u 2 = 3pT, J duP(u) (v-u) 2 = v 2 pT (18) 



we recast ( 17) into 



ir = 8v w +4v w 2 - 



t = ixa 2 pTt (19) 



Since V — v 2 , we have 

did d 2 

&V ~ 2v dv ' W 2 



Id Id 2 

4i) 3 dv 4v 2 dv 2 



Using these identities we re- write ( 19 1 as 



dr dv dv 2 
This kinetic equation admits the scaling solution 

1 



f(v,t) 



8-KT 3 



-v jr 



(20) 



(21) 



(22) 
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For instance, the average speed of the particle is 

/>oo 

(v) = / v f(v, t) Attv 2 dv — 3r 
Jo 



and more generally 



(23) 



The above analysis can be straightforwardly extended 
from three to d dimensions. The results up to ( 14 ) require 



obvious amendments, e.g. in equation (14) we must re- 
place a 2 by a d_1 . The integrals ( 15a H 15b) become (see 
Appendix [Bl 



J De(u-e)(g.e) = 4(u-g) 



Ve (u ■ e) 2 (g ■ e) 2 = ^— - (u ■ g) 2 + — 



where A, B are constants defined by integrals: 



B 



A 



De (g 



B 



De(g- 



(24a) 

g 2 u 2 
(24b) 



(25) 



The governing kinetic equation that generalizes Eq. (21) 
reads 



df 
dr 



■df , d 2 f 



dv 



dv 2 



t = 2a d ~ 1 ApTt. 



(26) 



Interestingly, in all dimensions the constant B drops from 
the final equation; the constant A is essentially irrelevant 
as it is absorbed into the new time variable r. 



Equation (26 1 is much simpler than Eq. ( 11 ) and it can 



be solved by employing the Laplace transform (see Ap- 
pendix [C]) . The asymptotic solution of Eq. ( 26 ) is again 
a pure exponential 



f = [Cl d T(d)Y 



-d e -v/r 



(27) 



is the area of the unit sphere in d 



'Vffl . . df df f f 

The constant in §27] has been chosen to = 4 I duF(u) g 1-7 / De (u • e)(g • e) 

JduP^g^J Be(u-e) 2 (g.e) 5 



where = 
dimensions. 

ensure the normalization: j dv/(v, t) = 1 

In two dimensions, Eqs. (p6|-|27]) have been derived 
in Ref. [10 in the realm of a stochastic model for Fermi's 
acceleration. Even earlier, the exponential velocity dis- 
tribution was found to occur in another stochastic model 
for Fermi's acceleration |19j in which a particle is bounc- 
ing in a container of fixed volume with boundaries de- 
forming in a chaotic manner. In this case, the velocity 
distribution becomes exponential independently of the 
container's shape and the deformation protocol. 



IV. MONOATOMIC GAS 



the particle and an atom separated by distance r can be 
described by a potential function U(r). In the long time 
limit when the particle velocity becomes large, only the 
small r behavior of the potential U (r) matters. In this 
limit, the repulsion part of the interaction dominates and 
it usually diverges algebraically in the small separation 
limit 



(28) 



as r — > 0. For example, A = 12 for the Lennard- Jones 
potential (in three dimensions). 

To estimate interaction size we can use the criterion 
U(r*) ~ q 2 , from which we find r* and the cross section 



area a 



n d-l. 



l/A 



(d-l)/X 



The term ga d ~ 1 De characterizing the hard-sphere gas 
should be replaced by the term ga*De in the general 
case. In one dimension, the interaction law is irrelevant 
and the problem reduces to the hard-core interaction. 
In higher dimensions, the Boltzmann equation depends 
on the interaction exponent A as it contains the factor 
.gcr* ~ g 1 ~ 7 with 7 = 2(d — 1)/A. In the long-time limit, 
the particle is very fast, so it is scattered only when it 
greatly approaches the atom, that is the separation is 
small and therefore the above analysis is asymptotically 
exact. Thus we must merely replace g by g 1-7 in the 
Lorentz-Boltzmann equation. This gives 



a/(y) 

dt 



duP^^-T/De [/( v ')-/(v) 



(29) 



where we absorbed the (rQe 1 l x ) d 1 factor into the time 
variable. 



To simplify the Boltzmann equation ( 29 ) we repeat the 



same steps as for the hard-sphere gas to yield 



a 2 / 

dV 2 



(30) 



where we have kept the terms with the first and second 
order derivatives in V as asymptotically they provide the 
leading contribution. Computing the angular integrals 
[as in Section III and Appendix fBI we arrive at 



1 df = d 2 f 
AA dt dV 2 



du P(u) v 



l-TT 2 2 



(u-v) 2 ] 



+ !£|duP(u) ff 1 - 7 (u.g) 



(31) 



Consider now a general case of a monoatomic gas. It 
is then natural to assume that the interaction between 



in the leading order. Thus the entire effect of the inte- 
gration measure is captured by one number, A. 
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To simplify the first integral on the right-hand side of 



(31 1 we write 



+ V- 1 -" 1 [(l- 7 )(vu) 2 + W 2 u 21 

where we have kept only the leading and the sub-leading 
terms. The integral over the leading term vanishes. Us- 
ing (18 1 and (20) we recast Eq. (311 into 



dr 



{d-l) 



Of , d 2 f 



di 



dv 2 



(32) 



where the modified time variable is given by [we addi- 
tionally put the factor (roe 1 /*)"* -1 back into the time 
variable] 



r = 2A{r Q e 1/x ) d - 1 pTt 



(33) 



Although one cannot [20] compute the factor A without 
knowing the integration measure, it is just a number that 
can be absorbed into the definition of the time variable to 



arrive at a universal kinetic equation ( 32 ) that depends 
only on the interaction exponent A. 

,1+7 



The form of equation ( 32 ) implies that r 



suggests a scaling ansatz 



/ = T- Ad <S>(w) 



A = (l+ 7 )- 



This 



(34) 



Plugging ( 34 ) into (32 1 we obtain an ordinary differential 



equation for <&(«;) which is solved to yield 



$H = C expj-A 2 ™ 1 ^} , C = 



n d T(Ad) 



(35) 



Thus the asymptotic growth, (v) ~ r , of the average 
speed and the scaled velocity distribution have univer- 
sal behaviors, the only parameters that matters are the 
interaction exponent A and the spatial dimensionality d. 

To exemplify the speed growth we note that in three 
dimensions 

r when A = oo (hard sphere gas) 
t 3 / 4 when A = 12 (Lennard- Jones gas) 
r 1 / 2 when A = 4 (Maxwell molecules) 

By definition, the Maxwell molecules (MM) interaction 
[2T] leads to the collision integral that is independent on 
the relative velocity. Equation (29) shows that the MM 
interaction is characterized by 7 = 1, so the interaction 
exponent is given by A = 2{d — 1). Interestingly, for 
the MM particle-atoms interaction, the average velocity 
experiences standard diffusion and the scaled particle ve- 
locity distribution is Gaussian. 

Let us now estimate the range of the validity of the 
above results if the particle mass m is small but finite: 
< m « 1. For a while, the evolution follows the zero- 
mass limit, but eventually the particle equilibrates with 



the background. The crossover to this regime occurs 
when the particle velocity becomes of the order of 

Y 

m 

In the earlier regime, t < t c , we have (v) ~ t 1 /^ 7 ). The 
crossover time t c is therefore estimated from 



(roe^Y-i^ 



1+3 



that is, 



tr 



(roe i/A ) 



i/A\-(<i-i), 



(36) 



The dependence of the crossover time t c on the gas den- 
sity and the mass of the particle is easy to appreciate. 
On the other hand, the dependence of the crossover time 
on the gas temperature is a bit surprising: 

1. When 7 < 1, that is A > 2(d — 1) implying that 
the potential is harder than the MM potential, 
the crossover time decreases as the temperature in- 
creases. 

2. When 7 > 1, that is A < 2(d - 1) implying 
that the potential is softer than the MM potential, 
the crossover time increases as the temperature in- 
creases. 



V. DISPLACEMENT OF THE IMPURITY 

We now turn to the spatial behavior of the impurity. 
We begin with a heuristic analysis. In one dimension, 
the mean- free path is I = the average speed grows 
as pTt [see Eq. ([6j|], and hence the time interval between 
collisions is At ~ p~ x /(pTt). This leads to an estimate 
for the total number of collisions during the time interval 
(0,0 



t 

At 



Ti 2 

I 2 



(37) 



The standard random walk argument tells us that a typ- 
ical displacement of the particle is given by 



Xtyp 



Tt 



(38) 



Hence the displacement exhibits a ballistic, x ~ t, rather 
than diffusive growth with time. Another unexpected 
feature of the growth law ( 38 1 is that the gas density p 



does not affect the asymptotic. 

The situation remains the same for an arbitrary di- 
mension d and an arbitrary interaction. Consider first 
the hard-sphere interaction. The mean-free path is £ ~ 
(pa d_1 ) _1 and the average speed is v ~ pa d ~ 1 Tt, see 
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Eq. ( 26 ) . Proceeding as in the one-dimensional case we 
find 



and therefore 



r typ 





df+ 


= 2pT 




d 2 f-' 


t Tt 2 


dt 


dv 


dv 2 


At ~ ~W 


df- 


= 2pT 


'df± 






dt 


dv 


dv 2 



and then proceed as in Sect. [II] to yield 

pv(U-f-) (43a) 
pv(U-f-) (43b) 



WW 



(39) 



The striking feature of this growth law is that the dis- 
placement is asymptotically independent on the density 
of atoms and their size. 

If the particle mass m is small but finite, < m <C 1, 



the growth law (39) holds up to the crossover time t c 



when the displacement becomes of the order of 



(40) 



while for t > t c the ballistic growth ( 39 1 switches to the 
diffusive growth 



(41) 



The above heuristic argument can be extended to the 
case when the particle-atoms interaction is described by 
a potential. At any time, the model is close to the hard- 
sphere case with effective radius of the order of . But 
since the displacement growth ( |39[ ) is independent on a 
in the hard-sphere case, it will be independent on at 
any given moment, and generally independent on the pa- 
rameters of the interaction potential ( 28 1 . Thus the dis- 
placement obeys the same growth law ( |39[ ) independently 
on A and d. 

We now turn from heuristics to exact analyses. To de- 
termine the second moment of the spatial distribution we 
first express it through the velocity correlation function 

(x 2 (t)) = f dti [ dta («(ti)«(t2)) 
Jo Jo 

= 2 f dt x f dt 2 (v(h)v(h)) (42) 

Jo Jti 

To evaluate (v(ti)v(t 2 )} let us consider the impurity par- 
ticle that starts at the origin with velocity equal to zero 
(initial conditions are actually irrelevant as we are in- 
terested in the long time behavior however this partic- 
ular choice makes the computation more compact). In 
this case the probability distribution for v\ = v(t\) is 
given by Eq. @. To deter mine the velocity distribution 
of v 2 — v(t 2 ) we must use v\ as the initial condition. 
The corresponding distribution function (i.e. the condi- 
tional probability) f(v 2 , ial^i, t\) satisfies a kinetic equa- 
tion which is different from ^ as the derivation of the 
latter assumes that the distribution function is symmet- 
ric, f(v) = f(—v). Generally we write 



/(«) = 



f+(v) v > 
/_(-«) v<0 



Subtracting (43b) from (43a) we see that the anti- 

(44) 



symmetric part 

#«) = /+(»)-/_(«) 
satisfies a closed equation 



90 
di 



-2pT 



d_0 
dv 



dv 2 



2pv(j> 



(45) 



(while for the symmetric part ip(v) — f+(v) + f-(v), we 
recover Eq. (JsJ) ) . The initial condition is 



(v,t = ti) = S(v - vi) 



(46) 



and the boundary condition, which follows immediately 
from the definition Eq. (44), is 



cf>(v =0,t)=Q 



(47) 



The initial-boundary value problem ( 45 )— ( 47 ) is non- 
trivial, yet in the interesting long time limit the governing 
equation (45) simplifies to ^ = —2pv(f> (since v >• VT), 
or equivalently j£ = —v<f>/T. Therefore 

4>(v, t\v 1} h) = 5{v - vi) e - v{ - T - Tl)/T (48) 

The velocity autocorrelation function can be presented 
in a rather compact form 



{VlV 2 ) 



dvi Vif(l) 
dvi wi/(l) 



— OO 

oo 



dv 2 v 2 f(2\l) 
dv 2 u 2 0(2|l) 



(49) 



Note that only the anti-symmetric part of /(2|1) con- 
tributes to the 2-points velocity correlation function. For 
the higher-points velocity correlation functions both the 
symmetric and anti-symmetric part appear alternatively. 
For example the 4-points velocity correlation function 
can be written as: 

(viv 2 v 3Vi ) = 2 |n jf dv t vA /(1)0(2|1)^(3|2)0(4|3) 

where ip( v j t\v 2 , t 2 ) satisfies Eq. (JsJ) with the symmetric 
initial condition ip(v, t — t 2 ) = 8(v — v 2 ) + 5{v + v 2 ). 

Substituting into ( |49| the results for /(l) = /(vi,ii) 
and 0(2)1) = (j>{v 2 ,t 2 WM) [Eqs. Q and {48} ] we get 



{V\V 2 ) 



dv\ — exp 

n 

2r 2 



1 



T 2 - Tl 

T 



[l + ^-rOn/Tp 
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Note that the equal times velocity autocorrelation func- 
tion (ti = t 2 = t) reduces to (v 2 (t)) = 2r 2 . This result 
directly follows from ^ thereby providing a useful check 
of the consistency of our calculation of the velocity au- 
tocorrelation function. Plugging the velocity autocorre- 
lation function into Eq. ( 42 ) we obtain 



fM*J dTlTl J Tl [1 + (r 2 - r x )r x /Tf 



(50) 



Computing the integral over r 2 yields 



1 

WTJo 



dn n<i- 



l 



[1 + (r 2 - n)n/T} 



The first integral J dri T\ provides the leading contribu- 
tion. Recalling that r = 2pTt we arrive at 



(x 2 ) ~ Tt 2 



(51) 



This asymptotically exact result confirms the heuristic 



prediction (38). 



One can also compute higher-order velocity correla- 
tion functions, e.g. (V1V2V3V4) , and use them to compute 
higher moments of the displacement. For instance, 

(a; 4 ) = 4! JJJJ dtidt 2 dt 3 dt4 {v^v-iV^) 
o<t 1 <t 2 <t 3 <t i <t 

These computations are very laborious, so we do not 
present them; we just mention that using this method 
we were able to compute the asymptotically exact fourth 
moment of the displacement, 



(x 4 )~5T 2 t 4 . 



(52) 



in one dimension. 

Finally we note that the above procedure can be gen- 
eralized to higher dimensions. Even in the case of the 
hard-sphere particle-atom interaction, however, the ex- 
plicit computations are quite unwieldy. 



VI. VELOCITY-POSITION DISTRIBUTION 



velocity distribution function /(v, t). In studying the dis- 
placement one would also like to use a governing equation 
for the density function N(r, t) as a starting point. Un- 
fortunately, there is no closed equation for the density 
function N(v, t). 

In the one-dimensional setting, the governing kinetic 
equation for the joint distribution f(x,v,t) reads 



v— -2 t(—+v — 
dt dx V dv dv 2 



df 



(54) 



The left-hand side of this equation is exact, yet Eq. (54) 



is already a simplified version of the Boltzmann equation 
as the collision term is only asymptotically exact, namely 
it is appropriate when v 3> y/T. As we mentioned ear- 
lier there is no closed equation for the density function, 
N(r, t). If one tries to integrate the kinetic equation (54) 
over v, the convective term leads to a current term, i.e. 
^ J dv vf(v,x,t) = J^J(x,t), so the density is coupled 



to the current. One can then deduce from (54) an equa- 



tion for the current, but it will involve the second moment 
J dv v 2 f(v, x, t). This procedure leads to an infinite hier- 
archy which seems intractable as (essentially) all infinite 
hierarchies. 



The kinetic equation ( 54 ) is a linear partial differential 



equation with two coefficients depending linearly on the 
velocity v. The most difficult term in Eq. (54), namely 



the convective term (vS/)f, can be further simplified in 
the long time limit when v 3> \PT. Indeed, since the par- 
ticle speed grows (on average) with a constant rate, the 
particle experiences numerous collisions during a time in- 
terval when its speed is almost constant. Then the prob- 
lem is akin to the standard Lorentz gas where the particle 
undergoes a simple diffusion. The separation between the 
time scale at which diffusion appears (few collisions) and 
the time scale at which the particle speed changes ap- 
preciably allows us to replace the convective term by the 
diffusion term of a standard Lorentz gas. In one dimen- 
sion, the diffusion coefficient is D = v/2p, see [5J. In the 
present case we can use the same formula. Thus Eq. ( 54 ) 
becomes 



df 
dt 



= 2pT 



df 

dv 



dv 2 



v d 2 f 
2~p dx 2 



The calculations of the moments of the displacement, 
e.g. the derivation of equation (52), through the veloc- 



ity correlation functions are very cumbersome. It seems 
hardly possible to succeed in deriving the next moment, 



61 T 3 i 6 , 



(53) 



relying on the velocity correlation functions. 

Therefore we employ different procedures that utilize 
a Boltzmann equation for the velocity-position distribu- 
tion /(r,v,i). This joint distribution function provides 
a complete description of the evolution of the impurity 
particle. Recall that in studying the velocity distribu- 
tion function we relied on a shorten description for the 



As usual, it is convenient to use r = IpTt as the time 
variable. Then the above equation becomes 



dj_ 
dr 



df 
dv 



d 2 f 



d 2 f 



dv 2 Ap 2 T dx 2 



(55) 



In Eq. (55) we tacitly assume that v > 0. This is ob- 



vious regarding the last term on the right-hand side as 
the diffusion coefficient must be positive (the correct ex- 
pression is D = \v\/2p). The form /„ + vf vv of the 
collision term also assumes (see Sect. [TTJ) that v > 0. 
There is no need to separately consider negative veloci- 
ties, it suffices to take into account the reflection symme- 
try f(x,v,t) = f(x,-v,t). 





FIG. 1. Shown are simulation results (see Sect. |IX| . Contour- 
plot (top panel) and 3D-plot (bottom panel) of the distribu- 
tion function F(X, V) for a gas of hard spheres in one dimen- 
sion. For any given velocity (position) the green (blue) dashed 
line shows the value of the position (velocity) for which the 
probability distribution has a maximum. 



In the long time limit, the joint distribution function 
f(x,v,t) should approach the scaling form 



f(x,v,t) 



1 



F(X,V), X 



X 



V = — (56) 



where = \/Tt and v* = t. The reflection symmetry 
with respect of the velocity and the displacement [55] 
allows us to limit ourself to the quadrant V > 0, X > 0. 



By inserting ( 56 ) into ( 55 ) we obtain 
dF OF d 2 F 



2F + X 



dF 



V 



V 



dX ' ' dV ' dV ' ' dV 2 
The normalization condition 



V 



d 2 F 
dX 2 



(57) 



dx 



can be re-written as 



dv f(x, v, t) = 1 



/>oo />oo 

/ dX dV F(X,V) = 1 
Jo Jo 



(58) 



This explains the factor 1/4 in the scaling ansatz (56) 



In higher dimensions, we limit ourselves to the case 
of the hard-core particle-atoms interaction. Then the 



governing kinetic equation reads 
df -^Ap T (d df 



df 
dt 



Or 



2a 



V dv 



d 2 / 
dv 2 



-DV 2 f (59) 



Equation ( 59 1 is again asymptotically exact in the large 



time limit when the typical particle velocity greatly ex- 
ceeds the thermal velocity, v \/T. In this limit, the col- 
lision term simplifies to the first term on the right-hand 
side of Eq. (59) and the convective term (v • V)/ can be 
replaced by the diffusion term — DV 2 / as the transport 
is asymptotically diffusion with velocity-dependent diffu- 
sion coefficient. More precisely, the diffusion coefficient 
is given by [55] 



D 



2dAa d ~ 1 p 



(60) 



with the amplitude A known in the case of the hard- 



core interaction, see ( B9 1 . Using again the modified time 
variable is r [which for hard-sphere particle-atom inter- 
action is given by r = 2Aa d ~ 1 pTt, see (26)], and taking 
into account the spatial isotropy we recast ( 59 1 into 



dr dv 



dv 2 



d{2Aa d - 1 p) 2 T 



dr 2 



ldf 
dr 



(61) 



A solution to Eq. approaches a scaling form 



f(r,v,t) = (tl d )- 2 (rVTf) d F(V,R) (62) 



with scaled spatial and velocity variables 



R 



Tt 



V 



(63) 



With the choice ( 62 ) of the scaling form, the normaliza- 
tion requirement 



Qar^dr 



n d v d - l dvf(r,v,t) = l 



becomes 



OO pOO 

dR I dVR 



d - 1 V d - 1 F(R,V) = l (64) 



Using (62)-(63) we transform (61) into 



2dF+RF R + VF V + dF v + VF VV 
] ^F R + F RR ) =0 



(65) 



This is a linear elliptic (recall that R > 0, V > 0) partial- 
differential equation. Despite of linearity, Eq. ( 65 ) is 



difficult since the coefficients in front of derivatives in 



Eq. ( 65 ) vary with V and R. 



We treat above equations by using different techniques. 
The standard technique relying on the Laplace and 
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Fourier transforms is the most powerful. In Sect. VIII 



we derive the major result for the scaled joint distribu- 
tion of the impurity particle in the hard-sphere gas: 



F(IL V) = ^ I <l$< 1 



(sin^)" <™ 



Further, the scaled density distribution reads 
,-iVds-K d d/2 n,j 



N{R) = C d / ds 



(coshs) d 
In particular, in one dimension 

1 



C,= 



(27r) d 



N{X) 



, R x = - X 
cosh Hi 2 



(67) 



(68) 



while in three dimensions the density is 

3s/3 (4i?| + 7T 2 ) tanh R 3 - 8R 3 nVZ 

N ^ = ~8 i? 3 cosh J R 3 ' Ra = ~2- R 

(69) 

First, however, we describe an approach based on the 
direct computing of the moments and guessing from them 
the spatial distribution. 



VII. MOMENTS 

The moment approach deals with the moments of the 
joint distribution rather than with the joint distribu- 
tion itself. The moment approach has been used in ki- 
netic theory throughout its history (see e.g. [T2J[T3]) as 
the governing equations are very complicated and sel- 
dom tractable. The moment approach has also been ap- 
plied [TU1 US] to the Fermi's acceleration mechanism. For 
instance, in Refs. |19j the authors computed the moments 
(v n ) for small n, guessed the answer [namely (23)] for an 
arbitrary n, showed that the guess is correct, and ob- 
served that the exponential velocity distribution has ex- 
actly the same moments. Generally if one succeeds in 
computing the moments, one still has to recover the dis- 
tribution that has such moments. This is not rigorous 
as at best we have infinitely many integer moments (or 
only even integer moments as in examples below) and 
we want to restore the entire distribution function. If 
the distribution function is analytic (the fact which is 
usually unknown, but believed to be correct), the distri- 
bution function can be uniquely determined by (infinitely 
many) integer moments, so restoring such function is a 
technical problem. 

Another problem is that since the number of moments 
is infinite, it is usually impossible to compute them all. 
Having computed a few moments one can try to guess the 
rest and to check the conjecture using computer-assisted 
exact calculations. We have succeeded in guessing all 
even moments of the spatial displacement in one and two 
dimensions, and in reading off the density in one dimen- 
sion. The moment approach is therefore not really sys- 
tematic and it involves a guess work. The strength of 



the moment approach is that one can easily compute the 
basic moments, e.g. even moments of the displacement 
(R 2 ) , (R 4 ) , (R 6 ) , etc., or mixed moments like (R 2 V 2 ), 
and arrive at important conclusions (like the existence 
of correlations between the velocity and the spatial dis- 
placement manifested by relation (R 2 V 2 ) ^ (R 2 )(V 2 )). 

In our problem we eventually derived more comprehen- 
sive results using standard techniques (see Sect. VIII[ ). 
Still, the moment approach has a future. Indeed it is 
more powerful nowadays than it ever was as the tedious 
calculations of the moments can be exactly performed by 
a computer and if the resulting moments admit a sim- 
ple expression through well-known sequences, there is a 
good chance to extract such an expression by using The 
On-Line Encyclopedia of Integer Sequences [25] . Since 
the moment approach is rarely used, we illustrate it here 
as in our situation where the moment approach clearly 
gives highly non-trivial results. We begin with the one- 
dimensional setting. 



A. One Dimension 

In this subsection we will present a very strong evi- 
dence in favor of the announced result ( 68 1 for the spa- 



tial distribution. To establish (68 1, we turn (57 1 into an 
infinite set of relations 

(i + j)Mij = fMij-i + i{i - l)Mi_ 2 , i+ i (70) 

for the moments 



poo poo 

M itj = I I dXdVX l V j F{X,V) 
Jo Jo 



(71) 



The relation ( 70 ) is valid for all i > 2, j > 0. 



Using (|70j) one can compute moments with small in- 
dexes; for instance, one can establish (51 )— ( 53 1 . Fig- 
ure ([2]) illustrates the structure of the quasi-recurrent 
equation (701 and the procedure to calculate the first 
few spatial moments. One finds that (X 2n ) = M2 n ,o can 
be expressed as a weighted sum of Mo,i, . . . , Mo, n . This 
sum is then computed using the identity 



Moj = {V 3 ) = 



dVi 



V>=j\ 



(72) 



We now demonstrate this in practice. Specializing ( 70 ) 
to — (2, 0) gives M2.0 = Mq,i = 1 which is identical 



to Eq. ([51]). Specializing (|70|) to = (2,1) yields 

3M 2jl = Af 2i0 + 2M , 2 (73) 
Taking then (i,j) — (4, 0) we obtain M4 = 3M 2 ,i, or 



.1/ 



4.0 



M, 



0,1 



2M 



0,2 



(74) 



which is equivalent to (52). Further, specializing (70) to 



= (6, 0), (4, 1), (2, 2) and using Q we obtain 

M 6)0 = 5M 0) i + 10M ,2 + 6M 0)3 - 61 (75) 
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FIG. 2. Schematic representation of how Eqs. ( 70 1 can be 
iteratively used to calculate all the moments Mij = {R l V J ) 
with i = even (red circles). The moments Mo,j are known for 
all j > 0. At the first step the known value of Mo,i allows 
us to calculate M2,o- At the second step the already known 
M2,o and Mo, 2 are used to calculate M21, see (731. At the 



third step we compute M^o through Ma,i. The moments Mij 
with i — odd (blue squares) cannot be calculated using this 
approach. 



which proves ( 53 1 . The fact that we have been able to 
reproduce the values of the spatial moments calculated 
using the velocity correlation functions (Eqs. ( 51 )— ( 53 )) 
supports the claim that the replacement of the convection 
term by the diffusion term in Eq. ( |54[ ) is asymptotically 
exact. 

The computed even moments (X 2n ) are all integers 
which look familiar; indeed, up to the sign they are the 
Euler's numbers 



(X 2n ) = (-l) n E; 



2n 



(76) 



The Euler's numbers E n appear in numerous combinato- 
rial problems, as well as in number theory, topology, etc. 
The Euler's numbers arc defined by the Taylor series 



1 



cosh(y) 



E n y n 

n>0 



(77) 



Note that all the odd-indexed Euler numbers are equal 
to zero, while the even-indexed Euler number have alter- 
nating signs. 

The evidence in the exactness of ( |76[) is overwhelm- 
ing — using Mathematica, we verified ( |76[ ) for all even 
moments up to (A" 1000 ). 



To establish (68 1 we start by extending the range of X 



to the whole axis and calculate the Fourier transform of 



N sym (X) = ±N(\X\): 

^sym(^) 



dXe 



-isX 



^sym (^0 



E 

n>0 
n>0 



(-l) n s 2n (X 2n ) 



(2n)! 



2n 77 

j2njT 



1 



coshs 



(78) 



where on the first step we have expanded e and taken 
into account that N sym (X) is an even function of X, while 



on the second and third steps we have used ( 76 ) and ( 77 1 , 
respectively. Since 



„-isX 



x dX cosh(7rX/2) 



1 



(79) 



we conclude that N sym (X) 
equivalent to Eq. (|68|). 



cosh s 

l/[2cosh(7rA72)] which is 



The moment relations (70) have helped us to deter- 
mine all even moments (X n ), yet they do not allow one 
to determine even the simplest odd moment (X). Us- 
ing the spatial density (68 1, however, we can compute 



this moment (more precisely it is equal to and it 

represents the dimensionless average displacement): 



(1*1) 



dX 



X 



i cosh(7rA72) 

where G is the Catalan constant 

1111 

12 ~ 32 + 52 - 72 + " 

Hence the average displacement is given by 



8G 

^2 



G 



0.915965594. 



8G 



Tt 



Similarly, one can compute an arbitrary odd moment 



\X 



2fc-l\ 



2 2fc+1 (2fc- 1)! 



T 2k 



E 

m>0 



(2m + 1) 



2 A- 



We can establish some qualitative and quantitative 
features of the joint distribution without having its an- 
alytical expression. For instance, if the joint distri- 
bution has allowed the factorization, that is if it had 
the form N(X)F(V), then the moments would satisfy 
(|A7 \V\ j ) = (|-X'| i )(|V| i ). This is not so, e.g. 

(X 2 V) _ 5 (X 2 V 2 ) _ 7 {X^V 2 ) _ 331 
(X 2 )(V) ~ 3 ' (X 2 )(V 2 ) ~ 3' (X^iV 2 ) ~ 75 

etc. Qualitatively, these results are not surprising — the 
larger separation from the starting position, the larger 
(on average) the speed of the particle is expected to be. 
Mathematically, this implies an inequality 



(|AT \v\ s ) 
(\X\i)(\V\i) 



> 1 



(80) 
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for all i, j > 0. This inequality is indeed obeyed in all in- 
stances where we were able to compute the moments, for 
instance when both indexes are sufficiently small. Using 
Eqs. ( 70 ) we have also computed a few infinite scries, e.g. 



(x 2n \v\) 

(X*»){\V\) 

(x 2 \v\ J ) 

(X*){\V\i) 

{x*\v\ j ) 

(X*){\V\i) 

(xW) 

m{\V\i) 



1 



\E- 



2n+2 I 



2n + l \E 2n \ 
2 

1 + -j 
3 J 



H H ¥ 

75 J 15 J 



> 1 



(81) 



794 



116 



40 



7 H 7 H 7" 

549 J 183 J 549 J 



Thus in these cases the inequality (80 1 is valid. 

The correlation between the velocity and the displace- 
ment of the particle shows that the knowledge of the 
velocity distribution F(V) and the density N(X) pro- 
vides a limited information about the characteristics of 
the particle — the joint distribution function F(X, V) 
is needed to provide a complete (in the realm of kinetic 
theory) description. 



B. Higher Dimensions 



The normalization condition ( 64 1 suggests to define the 
moments via 



dR / dVR 
i) Jo 



i-\-d— \xrj-\-d— 1 



V 3+a - x F(V,R) (82) 



Multiplying equation ([65]) by Ri+d-iyj+d-l and inte _ 
grating we arrive at the moment relations 

(i + i)M U3 = j(j +d- ljMjj--! 

i(i + d-2) 
+ - d 1 Mi- 2d+1 (83) 

We can now proceed as in the one-dimensional case. 



Namely using relations ( 83 1 , we can in principle exactly 
compute any moment (rt 2 ™) = M 2n ,o by expressing it 
as a weighted sum of Mo i, . . . , Mo n . Then we use the 
known expression for Mqj 



+d—l _ 



n.i+d) 

r(d) 



(84) 



which is computed with the help of Eq. (|27J). This pro- 
cedure gives 

(R 2 ) = d (85a) 

(i? 4 ) = (d + 2)(d+ |) (85b) 

(R 6 ) = d~ l {d + 2)(d + A){d 2 + 2d + f§) (85c) 

Using Mathematica, we have computed the moments 
(R 2n ) = M 2n , up to (R 1000 ) in two and three dimen- 
sions. A few of these even-indexed moments are listed 



n 


Id 


2d 


3d 


u 


1 


i 
i 


i 
i 


2 


1 


2 


3 


4 
6 
8 


5 

61 
1385 
50 521 


32 
3 
544 
5 

63 488 


55 
3 

1687 
9 

8651 


10 


35 

2 830 336 
63 


3 

5 047 691 
81 


12 
14 
16 


2 702 765 
199 360 981 
19 391512145 


357 892 096 


437 804 783 


231 

30 460116 992 
429 

26 862 763 900 928 
6435 


243 

16 325 727 605 
243 

6 868 768 364 827 
2187 



TABLE I. The moments (R n ) in one, two, and three dimen- 
sions for small even indexes. 



in Table I. In contrast to one-dimensional results (also 
presented in Table I), the moments are no longer integer; 
apparently 24] they are non-integer for all (even) n > 4. 

We tried to identify the sequence (R 2n ) — M 2n ,o with 
known sequences |25j . Since most known sequences are 
integer, one can seek M 2n ^ as a ratio of integer sequences. 
In three dimensions one can write (R 2n ) = M n /3 n . The 
sequence M n is integer, but it does not appear in [5S] . In 
two dimensions we were more lucky: Seeking M 2n $ as a 
ratio of integer sequences we arrived at 



ryZn+lj AT, 
(R 2 n) = i ^_ 



- 1 -!) 



1 



(2n) 



\B 2n 



+ 2 



(86) 



where Bk are the Bernoulli numbers 26]. The evidence 
in the exactness of ( 86 ) is overwhelming (we have checked 
it up to n = 500). 



C. Tail of the density distribution 

According to our definition of the scaled density dis- 
tribution N(R), it satisfies 



/>oo 

/ dRR d - 1 N(R) = l 
Jo 



(87) 



In one dimension, N = [cosh(7rA/2)] _1 , and therefore 
the tail of the distribution is 



N ~ 2e^ x/2 



when X — > oo 



This exact asymptotic leads to the conjecture that gen- 
erally in d dimensions the leading asymptotic is exponen- 
tial. More precisely, we assume that 



N ~CR C e-^ R when R 



(89) 



where we have augmented the controlling factor e~^ R by 
an algebraic pre-factor R c and the amplitude C. The 
parameters /i, c, C are dimensionless, so they can depend 
only on d. 

In principle, the moments 



(R 2n ) = dRR 
Jo 



2n+d-l 



N(R) 



(90) 
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depend on the entire density distribution N(R). In the 
n — > oo limit, however, the integral in Eq. (90) is chiefly 



gathered in the tail ol the distribution. Hence we can use 
the ansatz (89). Plugging it into (90) we get 



(R 2n ) 



C / dRR 2n+c+d ' 1 e-^ R 
Jo 
C 

„ t ~r~j T(2n + c + d) 

,,2n-\-c-\-d v ' 
A* 



(91) 



when n> 1. 



In two dimensions, Eq. (86) that yields even moments 



involves Bernoulli numbers whose asymptotic can be 
extracted from the celebrated Euler's formula relating 
Bernoulli's numbers with the values of the zeta function 
at positive even integers: 



B 2k \ = 



2 (2k) 



(27T) 



(92) 



Thus we recast (|86|) into 
(R 2n ) = 



., n . _ 2 3 " +1 (4" +1 - l) n!n! 2(2n + 2)! 

n + 1 (2n)! (27r) 2n + 2 



C(2n + 2) 



Using Stirling's formula, we simplify the ratio 

, 2n, 



n! n! 
(2n)! 



'2n\ 2n 



4wn 



'irn 
22* 



We also notice that ((2n + 2) 
fore asymptotically C(2ra + 2) 
moment (i? 2n ) approaches to 



(R 2 



1 ~ 2- 2n ~ 2 , and there- 
1 for n > 1. Thus the 



(93) 



irn (2n + 1)1 



^271+2 

in the n — > oo limit. On the other hand, in two dimen- 
sions the asymptotic prediction (91 ) based on the ansatz 

(94) 



( 89 1 can be re- written in the form 
(R 2n 



C 

2n+c+2 



(2n) c (2n + l)! 



where we used the well-known asymptotic [26] 
r(m + a) 



T(m) 



when m — > oo 



The asymptotics (|93|) and (94) would agree if 

C 



2 n+3 
^-2n+2 



2n+c+2 



(2n) 



We get (i = n/y/2 by matching the dominant exponential 
factors. Matching then the sub-leading algebraic factors 
we get c = 1/2. Matching finally the amplitudes yields 
C = 2 5 / 4 7r. Therefore in two dimensions 



N ~ 2 5/4 tt VR e'^/^ when R 



oo 



(95) 



3(4) 




FIG. 3. Plot of ^ 2 ^ 2 i+2^ ^ vs - n f° r the hard sphere gas in d = 



1,2,3. Using Eq. (971 we extract the controlling exponential 
factor e~^ dR of the density profile at large R and we confirm 
that Hd = \~K\fd in d = 1, 2, 3. 



The asymptotics in one and two dimensions make plau- 
sible that the controlling exponential factor in higher di- 
mensions is 



N - expj-^i?} 



(96) 



Thus N - e-" aR with ^ 3 = = 2.720699 in three 

dimensions. To extract /X3 we proceed as follows. Using 
Mathematica, we have determined the exact values of the 
moments (R 2n ) = M 2n ,o up to (i? 1000 ) in three dimen- 
sions. Hence we can compute the ratio of consecutive 
terms and compare the outcome with the prediction of 
Eq. (91). The latter becomes (in three dimensions) 



(R 2 



(M3) 2 



(R 2n + 2 ) {2n + c + 3) (2n + c + 4) 



(97) 



Thus the quantity (2n) 2 (R 2n ) / (R 2n+2 ) should converge 
for n -> 00 to (^ 3 ) 2 = 3tt 2 /4 = 7.402203. This is indeed 
in excellent agreement with our findings (Fig. [3). 



D. Correlations 

As in the one-dimensional case, both in two and three 
dimensions there are correlations between the position 
and the speed of the impurity particle. In this subsection, 
we present a few results for the three-dimensional case. 
One can compute (R^V 3 ) for even i and arbitrary j. For 
instance 



(R 2 V 2 



13 



(R 2 V 4 



17 



(R 4 V 2 



(R 2 )(V 2 ) 9 ' (R 2 )(V 4 ) 9 ' (R 4 )(V 2 ) 
etc. suggesting again that the inequality 



(R l V j ) 
(R l )(Vi 



> 1 



991 
495 



(98) 
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is valid for all i,j > 0. One can compute the left-hand 



side of Eq. (98) for arbitrary j and sufficiently small i: 
(R 2 V j ) 



W)W) 

{R & V ] ) 

WW) 



2 

1+ 9 J 



208 4 , 

7 H 7 

495 J 99 J 



27074 
45549 



3 + 



236 
2169 



r + 



40 
6507 



f 



(99a) 
(99b) 
(99c) 



For instance, let us establish (99a I. First, we specialize 



( 83 1 to d — 3 and i = 2 to yield 

{j + 2)M 2J = j(j + 2)Afa i3 -_i + 2A/ . J+ i (100) 



Using (84) an d set ting d = 3 we get Af 0j +i = |(i + 3)! 
and therefore ( 100 1 becomes 



M 2J =jM 3j -_i + (j + 3)0' + 1)! 



(101) 



The form of this recurrence suggests to seek M^j in the 
form M2.j = jl Nj. This transformation leads to 



A^-i + (j + 3)(j + l) 



(102) 



Solving recurrence ( 102 ) subject to the 'initial' condition 
A^o = 3 [this condition ensures that M 2i o = (R 2 ) = 3] we 
obtain 

1 i 
Nj = 3 + £(Z + 3)(Z + 1) = -(j + 1)0" + 2)(2j + 9) 



z=i 



Since (i? 2 VJ 

?2 



2,; 



j!7V, = i(j + 2 )!(2j + 9) and 



(R 2 ){Vi) = 3M 0tj = |(i + 2)!, we have 



i(j+2)!(2j + 9) 2 . 



811) in one dimension. 



thereby establishing (99a). Using similar reasoning we 
have derived (99b)-(99c), as well as analogous results 



The ratios (99a)-(99c) suggest that 

POO 

/ dRR 2+2i F(R,V) = e- y P t (V) 
Jo 



(103) 



with Pi(V) being a polynomial of V of degree i. We al- 
ready know that Po(V) = 1/2 in three dimensions. (Gen- 
erally P (V) = l/(d- 1)!.) Using |(99a])-([99c] we arrive 
at the following explicit results for the polynomials Pi(V) 
with i = 1,2, 3: 



Pi(V) 



1 

2 + 






17 


34 T T 


10 


18 H 


27 


457 


81 




162 


+ 



V 1 

266 

"81 



(104) 



V 1 



140 
243 



E. Monoatomic gas 

In the case when the particle-atom interaction has a 



power law tail ( 28 ) in the small separation limit, the joint 
distribution approaches a scaling form 

f(r,v,t) = (n d y 2 (r A Vf t)~ d F(V,R) (105) 



with scaled spatial and velocity variables 

v 



R = 



Tt 



V = 



-A 



(106) 



The analog of equation ( 65 ) reads 
(1 + K)dF+RF R + AVF V + U" 7 [(d - j)F v + VF VV ] 
+V \ F R + F RR ) = (107) 

Here D is a numerical factor which quantifies diffusion in 
the Lorentz gas where the particle-scatters interaction is 
given by ( 28 ) . 



Multiplying equation JT07l) by pj+d-iyj+d-i and in _ 



tegrating we arrive at the moment relations 

(* + Aj)Mi d = j(j + d-l- 7 )M lJ _ 1 _. 



i(i + d-2) 



(108) 



To the best of our knowledge, the value of the numerical 
constant T> is not known. 



VIII. JOINT DISTRIBUTION 



Here we derive the announced results (66)-(67) by em- 
ploying an approach based on the combination of the 
Laplace and Fourier transforms. It proves easier to deal 



with original kinetic equations (61) rather than with its 
scaled version. As a bi-product, we can also see that the 
solution approaches the scaling form. 

We begin again with the one-dimensional setting and 
show that the Laplace and Fourier transforms allow one 



to solve Eq. ( 55 ) for an arbitrary initial velocity distri- 



bution. Then we generalize to higher dimensions. 



One Dimension 



It is convenient to study Eq. ( 55 ) on the entire line 
— oo < x < oo while for the velocity will be taken posi- 
tive, < v < oo, as previously. Performing the Laplace 
transform in the v variable and the Fourier transform in 
the x variable, we find that the transformed joint distri- 
bution 



dxe lqx I dve~ vk f(x,v,r) (109) 

oo JO 
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satisfies 

! + (*'-« 2 >t = -^ V = m <110) 

This linear hyperbolic partial differential equation can be 
solved using the method of characteristics. The charac- 
teristics are the curves in the (k, r) plane which are found 
from 



dk 



= k 2 -Q 2 



Solving this differential equation we get 
k = -Qcoth[Q(£ + r)] 



(111) 



(112) 



where £ parameterizes different characteristics. Along 
a characteristics, that is keeping £ fixed, the governing 
equation (110) becomes 



dg 
<h 



^— const 



-kg 



(113) 



becomes 



Using ( |112| we express k via £ and r, so that Eq. (113) 
dg 



dr 



Qcoth[Q(£ + T)] 5 



whose solution reads 

g = sinh[Q(£ + r)] G(£) 



(114) 



(115) 



Specializing (112) and (115) to r = we get 



g (k, Q) = sinh(QC) G(£), k = -Q coth(QO 
so that 

5 o[-Qcoth(Q£),Q] 



G(0 



sinh(Q0 



(116) 



Combining ( 115 > ( [TT6] ) we arrive at the exact solution 
for the transformed joint distribution 



Using (112), we massage the ratio and rewrite the argu- 



ment of go to to transform (117) into 

k + Q tanh(s) 



g{k,q,T) 



1 



cosh s + j| sinh s ° I 1 + ^ tanh(s) 



where we have used the notation s = Qt which has been 



for any initial distribution 



used previously, e.g. in (79). This exact solution is valid 



9o(k,q) 



dxe lqx dve~ vk: f(x,v,r = 0) (119) 
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Consider now the simplest initial velocity distribution 



f(x,v,r = 0) = 6{x)5(v) 



(120) 



which corresponds to the initially stationary particle at 
the origin. The governing equation Eq. (55) is formally 
applicable if v 3> \ T (since the simplification of the col- 



lision integral in Eq. (JT) leading to Eq.(55) is valid only 
under this condition), but we are now more concerned 
with finding the simplest solution, in addition the ini- 
tial condition is asymptotically irrelevant. For the initial 
condition ( 120 ) we get go — 1 and the transformed joint 



distribution becomes 
1 



coshs 



^ sinhs 



s = Qr = qVft (121) 



The dependence on k in (121) is very simple, so we 



perform the inverse Laplace transform and obtain 



dq 
2n 



-iqx 



s e 



-Vs coth s 



• sinh i 



Note that the above formula already has the scaling 
form (for the initial condition (120) the scaling form es- 



tablishes instantaneously). Extending the variable V to 
the whole axis (this amounts to replace V —> \V\ and 
divide by 2) and re-writing the distribution in the mani- 

F(X,V) 



festly scaling form (f(x, v, r) 



see (56)) we get 



F(X,V) 



ds 



isX ■ 



- Vs coth s 



sinh i 



(122) 



Integrating in velocity, N(X) = f Q °° dV F(X, V), we ar- 
rive at the announced result (68). 



We could not compute the integral ( 122 1 in a closed 



form, so we determined it numerically. The results of the 
numerical integration (Fig. |4| are in excellent agreement 
with the results of direct simulations (Fig. [IJ . The excel- 
lent agreement between theory and simulations is further 
shown in Fig. [5[j6] and provides a verification of our sim- 
ulation scheme and shows that the replacement of the 
convection term by effective diffusion is indeed asymp- 
totically exact. 



B. Higher Dimensions 

The joint distribution f(r, v, r) is isotropic in r and 
v. It is convenient to explicitly assume the latter, so we 
want to find f(r,v,r). We define the Laplace-Fourier 
transform of this distribution through 



g(q,k,r)=n d J dre^ r J 



cW" 1 e~ vk f{r,v,T) 



(123) 

We limit ourselves to the hard-sphere interaction. Ap- 
plying the Laplace-Fourier transform to (61 ) we obtain 



dr + [k Q) dk 



-dk g 



(124) 
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FIG. 4. Contour-plot (top panel) and 3D-plot (bottom panel) 
of the joint distribution function F(X, V) for hard spheres gas 
in Id, Eq. (I22l. 




FIG. 5. Values of F(X, V) for hard spheres gas in Id along the 
lines of fixed V = 0.0035,1,2,3,4. The continuous lines are 
obtained from the numerical simulations (see Sec. IX I while 
the symbols represent the values obtained by computing the 
integral \122\. 



where we have used the short-hand notation 

-,2 



Q 



q- q 



d{2Aa d - 1 p) 2 T d(2Aa rf - 1 p) 2 T 




FIG. 6. Values of F(X, V) for hard spheres gas in Id along 
the lines of fixed X — 0.0035, 1, 2, 3. The data are obtained 
as explained in Fig. [5] 



instead of (113 )-( 114 1 we get 



dg 
(It 



= —dkg = dQ coth[Q(^ + r)]g 



^— const 

Integrating we find 

g = (smhm + T)}) d G(0 
while the general solution 

-d 



7(fc, q,r) = ( cosh s + — sinh s 



/ k + Q tanh(s) 
50 v l + ^tanh( S ), 



with s = Qt. For the simplest initial velocity distribution 
/(r,v,T = 0) = *(r)<J(v) (125) 
the general solution simplifies to 



g = ( cosh s + ~ sinh 



(126) 



As a check of this result we set q = 0. Then s = Qt = 
and limQ^o Q _1 sinh s — r, so that Eq. (1261 becomes 
g(k,q = 0,t) = (1 + rk)~ d which is exactl y th e Laplace 



transform of the velocity distribution [see ( C4 1 



Thus the joint distribution is the inverse Laplace- 
Fourier transform of (126). Performing the inverse 



Laplace transform of ( 126 ) in k is easy. Therefore the 



final answer is the inverse Fourier transform. Re-writing 
the result in the scaling form we arrive at the announced 
scaled joint distribution (66). Similarly we obtain (67). 



Equations (66) and (67) involve integrals of the kind 



J(R) = J dse- lV * s R <f>(s) (127) 



The characteristics curves in the (fc, r) plane are defined 
by the same equation (111) as in one dimension, while 



The integral J(R) is actually rotationally invariant, 
J(R) = J(R), which becomes clear by noting that 
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we can simultaneously rotate R and s. Using spheri- 
cal coordinates we write ds = £ld-i (sin 9) d ~ 2 s^ 1 ds d9 
where 9 is the angle between s and R (that is, we have 
s • R = sR cos 9). This allows us to reduce the d— fold 



with velocity u. At the instant of collision the particle's 
velocity changes. Thus the update rules are: 



integral (|127| to the double-fold integral 



j(R) = n d _ 1 / dss 4 - 1 ^) 

Jo Jo 



d9 (sin t 



-2 i\/d sRcost 



The integral in 9 is computable, so one actually reduces 
(1271 to a single integral. 



For example in two dimensions we have 
F{R,V) = 2 

and 



dss 3 J °(^ SR } ^-Vscoths 



/•CO 

N(R) = 2 1 ds 
Jo 



(sinh s) 2 

Jq(V2sR) 



(cosh s) 2 

while in three dimensions we obtain 



F(R,V) 



and 



ttR 



N(R) = 



ds g 4 sm(vggg) e _ ys coth s 
(sinh s) 3 

sin (V3sR) 



nR 



ds : 



(coshs) 



(128) 



(129) 



(130) 



(131) 



Computing the integral on the right-hand side of 
arrive at the announced result (69). The integra 



ing the joint distribution in d = 2,3 (Eq.(|128|)-(|l3"oj)) 



131 1 we 



s defin- 



were evaluated numerically and the resulting distribu- 
tions are qualitatively similar to the one shown in Fig. [4] 
for the Id case. 



IX. NUMERICAL SIMULATIONS 

In order to verify our theoretical results we have used 
different types of numerical simulations. 

The most straightforward numerical approach to check 
our theoretical results would be to perform a full molec- 
ular dynamics (MD) simulation. We are interested, how- 
ever, in the evolution of a single particle in a gas of back- 
ground atoms. The MD simulations are very inefficient 
to study such a situation since they keep track and up- 
date the positions and velocities of all the background 
atoms that are unnecessary to compute the quantities of 
interest. Whenever possible we turn to less costly com- 
putational method. 

For the hard sphere gas in one and two dimensions, the 
in-homogeneous Boltzmann equation was simulated by 
stochastically updating the velocity and positions of 10 6 
and 10 8 particles respectively. A particle with velocity 
v travels for a time At from the last collision covering 
a distance vAt before colliding with a background atom 



^n+l ~t~ At n 

r«+i = r„ + w n At n 

v„+i = v» + 2e[(u - v„) • e] 



(132a) 
(132b) 
(132c) 



Under the assumption already used in writing down 
the Lorentz-Boltzmann equation, the quantities At, u, e 
are random variables whose distributions need to be spec- 
ified in order to have a complete description of the tem- 
poral evolution. The velocity update rule ( 132c ) can be 



understood by analyzing the collisions in the reference 
frame of the background atom (which in our case coin- 
cides with the center of mass reference frame). The key 
feature of the hard-sphere interaction is that the collision 
rate is proportional to the absolute value of the relative 
velocity g, so that the particle more often collides with 
atoms moving in direction opposite to its own. 

The random variable At is the first collision time which 
is distributed according to a Poisson process. This can be 
understood in the following way. The particle can collide 
with any background atom. The probability that the 
particle has not collided with the background atom i th 
up to time t is called Si(t). The survival probability Si(t) 
is decaying in time and satisfies a very simple differential 
equation: 



dSjjt) 
dt 



Si(t), 



(133) 



The rate of collision, r$, is proportional to the absolute 
value of the relative velocity with respect the i th atom. 
The probability that the particle has not collided with 
any atom up to time t is S(t) = YliLi &i{t)i where N is 
the total number of background atoms. Using Eq. ( 133 ) 
and the definition of S(t) we obtain 



dS(t) 
dt 



-rS(t), 



N 



(134) 



whose solution is a simple exponential decay with rate 
r. Note that S(t) is also the probability that the first 
collision happens at time t, i.e. S(t) is the distribution of 
the first collision time. Reintroducing the dependence on 
the particle velocity explicitly we obtain the probability 
P(At\v) that the particle with velocity v collides for the 
first time at time At: 



P(At\v) = r(v) exp(-r(v)Ai) 



N 



"(v) = ]Tr s ;(v) = 2ap(|v-u|) u 



(135a) 
(135b) 



Here ((-))u denotes the average over the velocity distri- 
bution of the background atoms, a is the radius of the 
hard-spheres and p is the number density of background 
atoms. The last equality in ( 135b ) has been specified for 
the two-dimensional case. 



17 



The probability of making the first collision with the 



r-j(v) _ |v - Uj| 
r(v) ~ iV(|v-u|) u 



(136) 



This equation can be understood in the following way. If 
it was equally likely to collide with any atom only the 
factor 1/N would appear in Eq. (136). The correction 
( ( v-ul) ) m (1^6) to this simple behavior describes 



the fact the the particle collides preferentially with atoms 
moving in direction opposite to its own. It is worth noting 
that this correction approaches 1 if v ^> {u) . 

The calculation of the total rate is difficult in any di- 
mension d > 1. It can be approximated by 



v-u 



u 



T 



(137) 



Only the limiting behavior for u > vT and v <C \T 
of Eq. (137) are important. We are interested in the 
large time limit when v 3> (u) and r(v) ~ v| . Equa- 
tion (137) correctly reproduces this limit. Moreover, 



Eq. (137) ensures that a particle with an unexpected low 
velocity (in the extreme case v = 0) will collide with a 
background atom with a rate proportional to the thermal 
velocity o f the b ackgr ound gas. 

Using ( 135a )-( 137) one computes the collision time 



At. Then a background velocity u is generated from the 
Maxwell-Boltzmann distribution ^ and it is accepted 
with probability /jj^z^T) (see Eq. ( |136| )). Finally the ran- 
dom variable e is generated from the distribution (13). 



The velocity distribution is in excellent agreement with 
the exponential scaling form. The density profiles are 
shown in Fig. [7] In one dimension, there is a perfect 
agreement with the theoretical prediction, Eq. (68). In 



two dimensions, the numerical simulation correctly re- 
produces the known values for the moments (R 2n ) (see 
Table |l| and agrees with the prediction (196J) for the tail. 

In the one-dimensional case, every velocity distribu- 
tion of the background atoms is stationary (since in a 
two-body collision the atoms merely exchange their ve- 
locities). In particular it is possible to chose a uniform 
velocity distributio n for — u max < u < u max . In this case 
the total rate (Eq. (135b|) can be calculated exactly and 
Eq. (136) can be enforced very efficiently. In this situ- 



ation we were able to stochastically update the velocity 
and positions of 10 s particles which allowed us to sim- 
ulate the joint distribution F(X, V) (see Fig. [IJ. It is 
interesting to note how the exponential character of the 
speed distribution F(V) is also present for F(V, X = 0). 
In the same way the character of the density distribu- 
tion N(X) persists for F(V = 0,X). In the contour-plot 
(top panel of Fig. [T]) we observe that the equiprobability 
line always cross the V-axis perpendicularly while they 
cross the X-axis at acute (obtuse) angle for X < X c 
(X > X c ) where X c ~ 0.8. This has the consequence 
that for any given velocity the maximum probability is 
always at X = (green line in Fig. [IJ while for fixed X 



the maximum probability is at V = only for X < X c 
(blue line in Fig. [IJ . The numerical result clearly show 
the lack of factorization: The joint distribution F(X, V) 
is not a product of functions of X and V. 

In two dimensions, we have also used a "brute-force" 
molecular dynamics simulations to investigate the case 
when the atoms interact between themselves and with 
the particle through the potential U ~ r~ A . This sim- 
ulation schemes is much more time-consuming than the 
stochastic update of the position and velocity of the par- 
ticle. For this reason we were able to simulate only 10 4 
particles. This is sufficient to check the scaling of the 
average velocity and displacement with time, but does 
not allow us to check the full distribution. In our system 
the background atoms arc affected by other atoms and 
insensitive to the presence of the particle; the particle is 
affected by the atoms. Computationally this property is 
implemented in a simple way. At each time step of the 
molecular dynamic simulation we calculate the total force 
acting on a background atom summing only the contri- 
butions from the other atoms (no contribution from the 
particle). The total force acting on the particle is ob- 
tained summing all the contributions from the atoms. 

Numerically it is convenient to simulate many indepen- 
dent particles in the same background gas of atoms. Usu- 
ally, even if the particle-particle potential is set to zero, 
particles interact indirectly via the background gas. In 
our case, the particles do not affect the background atoms 
and are totally independent from each other. We have 
simulated 10 independent particles in the same back- 
ground gas of 5 • 10 3 atoms. For the reason explained 
before this simulation scheme is equivalent to 10 4 runs of 
a single particle in a background gas of 5 • 10 3 atoms. 

The equations of motion have been numerically inte- 
grated using the velocity- Verlet algorithm [28] . The time- 
step of the numerical integration was reduced during the 
time evolution in order to keep the average particle's dis- 
placement during a single time step constant and smaller 
than the mean free-path of the gas. The initial positions 
of the background atoms and of the particles were ran- 
domly drawn from the uniform distribution inside the 
simulation box with periodic boundary conditions. The 
initial velocity of the particles were drawn from the dis- 
tribution S(v — vo)/2n while the initial velocity of the 
atoms were generated from the Maxwell-Boltzmann dis- 
tribution and were rescaled in order to ensure that the 
total energy (~ T) of the background gas had a fixed 
value. 

The results of different simulations at fixed density and 
fixed interaction exponent are shown in Fig. [8] and [9] re- 
spectively and are in excellent agreement with theoretical 
predictions. 



Finally, the quasi-recurrent relation ( 108 1 has been it- 
eratively solved (as shown in Fig. [2] and explained in the 
text) with Mathematica to calculate exactly the moments 
of the spatial distribution (R 2n ) up to n — 500 for the 
hard sphere gas in d = 1, 2, 3. 

In Fig. [3] we show the ratio (2n) 2 (R 2n ) / (R 2( - n+1 1) 
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FIG. 7. Density profile for the hard sphere gas vs. the rescaled 
variable R = r/yTt. The numerical simulations in d — 1,2 
are compared with the theoretical predictions, Eq. (68 1 and 



Eq. (129 1 (inte grat ed numerically). The theoretical prediction 
for d = 3, Eq. (69 1, is also shown. 



<r> for h 
<r> for ^=6 
<r> for h=4 
<r> for X=2 
<v> for A.=8 
<v> for X=6 
<v> for h=A 
<v> for X=2 




FIG. 8. Average particle velocity and displacement in two di- 
mensions with a particle-atom interaction potential diverging 
as U ~ r - \ In all cases the density of the background gas is 
p = 25%. The slopes of the fitting curves (dashed lines) are 
0.5, 0.66, 0.74, 0.79 (bottom to top), all in excellent agreement 
with the theoretical prediction A/(A + 2). The solid line has 
slope 1 and it is a guide for the eye. 



which allows us to extract the asymptotic exponential 
decay of the density distribution. 



X. SUMMARY 

We have analyzed the behavior of a very light particle 
in an equilibrium background gas. We have shown that 
in the long-time limit, the average particle displacement 
grows linearly with time and proportionally to the ther- 
mal velocity of the background atoms — the density of 
the gas, the size of atoms, and the details of the interac- 
tion between the particle and the atoms do not affect the 
asymptotic. The average particle velocity also grows in 
a rather universal way and the scaled velocity distribu- 




<r>forp=10% 
<r> for p=25% 
<r> for p=40% 
<v> for p=10% 
<v> for p=25% 
<v> for p=40% 
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time 



FIG. 9. The same system as in Fig. |8| with fixed interaction 
exponent A = 8 and varying density. The slope of the fitting 
curves (dashed lines) is A = 0.79 in all cases, while the in- 
tercepts are b — 0.79, 1.60,2.41 (bottom to top). Note these 
values are in the ratio 1 : 2.02 : 3.05 in excellent agreement 
with the theoretical prediction ( |26[ ) bi/bj = {pi/pj) A which 
gives 1 : 2.08 : 3.03. 
guide for the eye. 



The solid line has slope 1 and it is a 



tion approaches a scaling form which is generically non- 
Gaussian (the only exception is when the particle-atoms 
interaction is described by a Maxwell potential). 

For the hard-sphere particle-atom interaction in ar- 
bitrary dimensions, we have computed the asymptot- 
ically exact velocity distribution, position distribution 
and joint velocity-position distribution. The most com- 
plete results for the joint distribution have been derived 
using a combination of Fourier and Laplace transforms. 

In one dimension, we have also determined the prob- 
ability density for the particle displacement using a 
less standard moment approach. Specifically, we have 
guessed an exact expression for the moments (r 2n ), which 
we verified by exact (Mat/iematica-assisted) calculations 



of the moments up to 



,.1000 



), and we found the prob- 



ability density that results in these moments. We have 
also guessed an exact expression for the moments (r 2n ) 
in two dimensions and we have confirmed to the same 
depth as in one dimension. Further, we have used the 
moments to establish the large displacement tail of the 
probability density and to study the correlations between 
the velocity and displacement of the particle. 

Our theoretical predictions are in perfect agreement 
with the numerical simulations providing strong evidence 
that our simulation scheme is correct and that the sim- 
plification of the collision integral and the replacement 
of the convective term by effective diffusion are indeed 
asymptotically exact in the limit when the particle ve- 
locity greatly exceeds the thermal velocity of atoms. 

The Lorentz model was originally suggested [T] as an 
idealized model of electron transport. Quantum mechan- 
ics is of course essential for this problem. In the context of 
the quantum particle in a container of fixed volume with 
boundaries deforming in a chaotic manner (a stochas- 
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tic model for Fermi's acceleration of the quantum par- 
ticle), some mostly numerical work has been done (see 
e.g. Perhaps the most interesting extension of the 

present work is to analyze the quantum version of our 
model. 
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Appendix A: Approach to Scaling 

In one dimension, atoms merely exchange their veloc- 
ities, so there is no relaxation and any velocity distribu- 
tion P(u) can be taken as an equilibrium distribution. As 
an example, consider the bimodal velocity distribution 



P(u) 











- 




+ s( wH 


-I) 











(Al) 



(The bimodal distribution is often used in studies of 
the one-dimensional Boltzmann equation, see e.g. [16].) 
Note that for the bimodal velocity distribution the condi- 
tion of Eq. ([3]) holds; further, the density and the temper- 
ature of the background gas are p = 2,T = 1/4. There- 
fore r = 2pTt — t and the scaling solution Q becomes 



f(v,t) 



-\v\/t 



2t 



Let us now try to establish exact results starting with 
initial condition 



f(v,t = 0) = S(v) 



(A2) 



The velocity distribution cannot approach the smooth 
distr ibution Q. For the bim odal velocity distribution 
(Al) and the initial condition (A2), the particle velocity 



can be only integer: 

00 

/(«,*)= ]T Pn(t)S(v- 



(A3) 



The amplitudes P n (t) are still expected to behave as 

1 



Pn{t) 



-\n\/t 



2t 



(A4) 



in the limit \n\ — > 00 and t — > 00, with n/t being finite. 

To probe the exact behavior we insert ( Al I and ( A3 ) 
into the Boltzmann equation ([T]) and deduce an infinite 
set of rate equations 

Pn = ( n - i J P n -i + (n + I) P n+1 - 2nP n (A5) 



for n > 1 and 

P = Pi - Po (A6) 
(It suffices to consider P n with n > 0; with initial condi- 



tion ( A2 ) , the particle velocity is a manifestly even func- 



tion of v and therefore P- n = P n . 



To treat (A5|- (A6) let us use the generating function 

(A7) 



Note that 



nt,z) = P (t) + 2j2Pn(t) z n 



(A8) 



explaining why we have chosen the definition ( A7 1 of the 
generating function instead of J2 n>0 P n (t) z n . 



Utilizing the generating function ( A7 ) we recast an in- 
finite set of rate equations ( A5 1- ( A6 1 into a single partial 
differential equation 



~dt 



d 1 ? ^ (1 - zf 
dz 



2z 



Po(t) 



1 



2z 



9 (A9) 



We want to solve (|A9| subject to the initial condition 

(A10) 



Pn{t = 0) = <5„.o, or equivalently 

y(t = 0, z) = 1 

and the boundary condition (lASh . 



Using C= 1/(1 — z) instead of z, we re- write (A9) as 
Po(t) 



d3> 

at 



1 - 2C 



d( 2C(C - 1) 2C(C-i) 



IP 



(All) 



The t ransformation £ = (t + Q)/2, t] = (t — ()/2 recasts 
(lAllh into 



dr) 



f , o(C + 7 ? ) + [i-2(g-7 ? )]y 

2(C-»7)(e-»7-l) 



(A12) 



To solve ( A12 ) we note that its homogeneous version, 

dy = 1 - 2(e - /?) y 

dv 2(e-r?)(e-r?-l) ' 
has a general solution 

where Q(l;) is an arbitrary function of ^. Then a solution 
to the full equation (A12) can be sought using the varia- 



tion of constant technique. In the present case we must 
actually vary the function Q(£), namely, we should seek 
a solution of the form 



nt,v) = \/(S-v)(Z-v-i)Q(Z,v) 



(A13) 
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Plugging ( |A13 ) into (A12) we obtain a simple equation 
for Q which is integrated to find a final solution. Return- 
ing back to the variables (t, Q we get 



nt,o = vc(c-i)Q(*,o 



with 



Q = 



1 



V(* + 0(* + C-i) 

1 ^ . Po(r) 



(A14) 



(A15) 



2 ./o [(* 



+ C)(<-t + C- 1)] 3/2 



Equations (A14)-(A15) give rather formal results as 



we haven't yet extracted Po(t). However, on this stage 
we can already confirm the emergence of scaling (A4|. 



Indeed, assuming that Po{t) decays and approaches to 
zero as t — ¥ oo, we conclude that the integral term on 
the right-hand side of (A15) is asymptotically negligible 
and therefore Q ~ l/(t + Q. Therefore (A14| becomes 
y ~ Q/(t + £), where we additionally consider the large C, 
limit. Hence T ~ 1/(1 + i/C) = l/(l + t-tz). Expanding 
this result we get 



Pn(t) = 



1 t n - Y 

2(1 + *)' 



which in the scaling limit n — >■ oo and t — > oo with n/t 
being finite is indeed equivalent to (A4). 



Appendix B: Angular Integrals 



Let us first prove the validity of relation ( 24a I with A 



defined in (25). The integral in (24a) is equal to (J • u), 
where J = jT>e (g ■ e) e. Due to symmetry, the vector J 
must be directed along g. Hence 



Ag 



(Bl) 



where the amplitude A is independent on g since J scales 
linearly with g. Computing the scalar product of g and 
J we obtain 



A=l(J-g) = l f 2)e(g-e) 2 

g 9 J 



(B2) 



Using (Bl I we arrive at 



J De(u-e)(g-e) = (u- J) = A(u ■ g) 



which together with (B2) lead to (24a 



To establish ( 24b I with B defined in ( 25 ) we note that 



the integral in Eq. ( 24b ) is equal to (u • T • u) , where 

T = J De(g-e) 2 ee (B3) 

Tensor T depends only on vector g, so it must read 

T = C 1 gg + C 2 g 2 U (B4) 



where U is the unit tensor. To determine the amplitudes 
C\ and C2, we compute the trace of tensor T and the 



product (g • T • g). Using (B4) we find 

Tr(T) = (d + dC 2 )g 2 
(g • T • g) = (Ci + C 2 )g 4 



If instead we use ( B3 ) we get 



Tr(T) = J De(g-e) 2 = Ag 2 
(g-T-g)= J De(g-e) 4 = Bg 4 



(B5a) 
(B5b) 



(B6a) 
(B6b) 



where we have used the definitions of A and B, see (25 1 



and C2 via A and B: 



Comparing ( B5 ) with ( B6 ) we express the amplitudes C 



Ci 



dB — A 

d-l 



C 2 



A-B 



(B7) 



yielding indeed ( |24b ). 

For the three-dimensional hard-sphere gas, the integra- 
tion measure is given by Eq. ( 13 ) and therefore 



A = 1 j d 2 e 0( S ■ e) (g • e) 3 
B=^Jd 2 e 6(g • e) (g • ef 



(B8) 



Let us now introduce spherical coordinates with the axis 
along g. We have d 2 e = 2ir sin $ eft?, (g • e) = g cos$; 
the term #(g • e) limits the integration over the range 
< ■& < tt/2. Thus 



A = 2tt 



it/2 



sin 



and similarly B = 7r/3. Thus we obtain ( |15a )-( 15b I. 
(See Ref. [3] for the computation of integrals similar to 
(15); such integrals often appear in kinetic theory of the 



hard-sphere gas.) 

For the d— dimensional hard-sphere gas, we have the 
same expression (B8) for A and B, the only difference is 
that de = fld-i (sini9) d ~ 2 d-d. Computing A yields 



A = n d _i 



*V2 J n (d-l)/2 

( s ind) d ~ 2 (co S d) 3 dV= T ^ m ^ (B9) 



Appendix C: Exact Solution of Eq . (|26[ ) and Analysis 
of Solutions of Eq. |32|) 



Let us first solve Eq. (26) using the Laplace trans- 
form. Note that in Eq. (26) the variable v varies in the 



range (0, +00) and therefore we use the Laplace trans- 
form rather than e.g. the Fourier transform. In any 
number of dimension we define 



g(k,r) = fld 

Jo 



dv v 



d—l „—vk 



f(v,r) (CI) 
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where £l d = r(d/2) 1S ^ e area °f the unit sphere in d 
dimension. According to this definition, the function g 
satisfies the boundary condition g(k = 0, r) = 1 and the 
initial condition 



g (k) = g(k,T = Q) = n d I dvv d - 1 e- vk f(v,T = 0) 



Applying the Laplace transform to Eq. ( 26 ) yields 

-dk g — k 2 



dg 

dr 



, 2 dg 
dk 



(C2) 



The right hand side can be rewritten as — k 2 d -^(k d g) 
thereby suggesting to use the function h = k d g instead 
of g. One gets h T = —k 2 hk, or equivalently 



dh _ dh 
dr dn ' 



k = k 



(C3) 



A general solution to the simple wave equation (C3) is 



h(K, t) — H(k + t) where H is determined by the initial 
condition: H(k } t = 0) = H{n). Returning to the orig- 
inal function g we arrive at the general solution for the 
Laplace transform 



](k,r) = (l + Tk)- d g a 



1 + kr 



(C4) 



As an example of the initial distribution with a com- 
pact support (that is, vanishing for sufficiently large ve- 
locities) consider the isotropic distribution with fixed ini- 
tial speed vq. In other words, let 



f(v,T = 0) = 



S(v - V ) 



(C5) 



In this case [3U] the solution reads 

1 



(l+Tk) d 



exp 



v k 



1 



Expanding the exponential and separately performing 
the inverse Laplace transform of each term we obtain 



f{v,r) 



1 



i-vo/r) n 



n d r(d)T d ^ 



1 F 1 [n + d-d-- v -} 



where \F\ is the confluent hypergeometric function. The 
asymptotic behavior (r 3> vq) of (C6) is given by the 



first term (n = 0) in the sum and is equal to 



f(v,r) 



n d T(d) r d 

where we have used the identity iF\[d\ d; z] = e z 



(C7) 



As an example of an initial distribution with infinite 
support, consider an exponential distribution 



f(v,r = Q) = 



1 



-v/vq 



(C8) 



n d r(d) v d 

In this case, the velocity distribution remains exponential 
throughout the evolution 



f(v,r) 



1 



-v/(v +t) 



n d T(d) (vq 



(C9) 



The asymptotic (r 3> z;o) behavior of the solution (C9) 
is again given by (C7). 



These two examples illustrate the general behavior 
which can be deduced from the general solution (C4|: 



If the initial velocity distribution decays exponentially or 
faster, the asymptotic behavior of the velocity distribu- 
tion is universal (that is, independent on the initial veloc- 
ity distribution) and given by ( |C7[ ). If the initial veloc- 
ity distribution decays slower than exponentially in the 
v — > oo limit, the long time asymptotic behavior is given 
by Eq. (C7) apart from the tail region. For instance, if 
f(v,r = 0) ~ v~ v as v — > oo, the asymptotic velocity 
distribution is given by (C7) when < v (v — d)r In r, 
while for v 3> (y — d)r In r the initial distribution domi- 
nates: f(v,r) ~ v~ v . 

Essentially the same qualitative behavior is valid in 
the general case of the potential particle-atom interaction 
(28). The governing kinetic equation (32) describing the 



long time behavior is substantially more difficult than 
Eq. ( 26 ) corresponding to the hard-sphere interaction, 



e.g. applying the Laplace transform to Eq. (32) does not 



lead to a closed equation for g(k, t). Therefore it is much 
harder to prove rigorously that the asymptotic is given 
by (34)-(35). A non-rigorous, but physically convinc- 



ing, argument relies on the existence of a one-parameter 
family of exact solutions generalizing the scaling solution 
([34])— ( 35 ) . Indeed, let us start with an initial velocity 
distribution l30l 



f(v,T = 0) 



C 
7/i 



exp< 



-A 



Vq 



I/A" 



(CIO) 



where vo is a parameter and the constants C and A are 
the same as in Eqs. (|34|-([35|. A solution of Eq. ([32l 



subject to the initial condition (CIO) reads [3T] 



f = C 



1 1 a \ —Ad 

i/M „ ) A 2 



exp< -h 1 



Obviously, the velocity distribution (Cll) approaches 
the scaling form (34 ) — ([35]) in the long time limit. This 




(Cll) 



strongly suggests that for an arbitrary initial velocity dis- 
tribution that decays as exp{— const, x w 1 ^} or faster, 
the asymptotic behavior is given by p4|)~(|35"|). For 
the initial velocity distribution decaying slower than the 
above stretched exponential, the asymptotic velocity dis- 
tribution is still given by Eqs. ( 34 1 — (35) in the major 



range and only the tail region is dominated by the initial 
velocity distribution. 



